% function test_possion_p
clear all;


caseNum = 4; epsilon = 0.005;

disp_d = 3; tri_temp = template_mesh_tri(disp_d);

[p,e,t] = initmesh('lshapeg','Hmax',0.35,'Hgrad',1.99);
V = [p(1,:)' -p(2,:)']; T = t(1:3,:)';

max_err = 5e-3;   % 5e-3 5e-4 5e-5 
[T, E, ET, TE] = build_fem_mesh(V, T);
[V, T, E, ET, TE] = refine_L_shape(V, T, E, ET, TE, 12);

d = 3*ones(size(T,1),1);
d_max = 15;  
[cr,desc,asce] = pattern_matrix_BB(d_max);

Int_UU_cell = cell(d_max,1);
Int_UV_cell = cell(d_max,1);
Int_VV_cell = cell(d_max,1);
Int_cell = cell(d_max,1);
% data for posteriori estimator
    Int_UU2_cell = cell(d_max,1);
    Int_UV2_cell = cell(d_max,1);
    Int_UW2_cell = cell(d_max,1);
    Int_VV2_cell = cell(d_max,1);
    Int_VW2_cell = cell(d_max,1);
    Int_WW2_cell = cell(d_max,1);
for degree = 1:d_max
    Mat = build_dim2(degree-1);
    Mat2 = build_dim2(degree-2);
    Id = diag(ones((degree+1)*(degree+2)/2,1));
    Du = degree*de_cast_step(Id,degree,1,0,-1,desc); % the direction derivate on v2-v1
    Dv = degree*de_cast_step(Id,degree,0,1,-1,desc);  % the direction derivate on v3-v1
        Duu = (degree-1)*de_cast_step(Du,degree-1,1,0,-1,desc); % the direction derivate on v1-v3
        Duv = (degree-1)*de_cast_step(Du,degree-1,0,1,-1,desc);  % blending partial derivatives
        Dvv = (degree-1)*de_cast_step(Dv,degree-1,0,1,-1,desc); % the direction derivate on v2-v3
        Int_UU2_cell{degree} = Duu'*Mat2*Duu;
        Int_UV2_cell{degree} = Duu'*Mat2*Dvv;
        Int_UW2_cell{degree} = Duu'*Mat2*Duv;
        Int_VV2_cell{degree} = Dvv'*Mat2*Dvv;
        Int_VW2_cell{degree} = Dvv'*Mat2*Duv;
        Int_WW2_cell{degree} = Duv'*Mat2*Duv;
    Int_UU_cell{degree} = Du'*Mat*Du;
    Int_UV_cell{degree} = Du'*Mat*Dv;
    Int_VV_cell{degree} = Dv'*Mat*Dv;
    Int_cell{degree} = build_dim2(degree);
end


flag = 1; count = 1; tic;
while flag    
    flag = 0; % solved again , then keep the flag 0  
    fprintf('\n\nStep %d : Sorting dofs... ',count);
    t0=cputime;  % sort the global dofs initially.
    [dofs,n_dof] = distribute_dof_spline(d); % init and reserve spaces for dofs and n_dofs
    t1=cputime;
    fprintf('\nFinished! Spend %3.20f secon',t1-t0);
    
    fprintf('\nAssembling %d elements...',size(T,1));
    t0=cputime;
    
    M = mass_bb_2(dofs, V, T, d, Int_cell);
    A = stiff_bb_2(dofs, V, T, d, Int_UU_cell,Int_UV_cell,Int_VV_cell);
    F = bnet_2(dofs,V,T,d,'elliptic_f',caseNum,epsilon);
    b = M*F;
    [B,bf]= elliptic_bc(dofs, V, T, TE, ET, d,'elliptic_g',cr,caseNum,epsilon);
    H = cr_matrix(dofs,V,T,TE,ET,desc,asce,cr,d,0);
    
    t1=cputime;
    fprintf('\nFinished! Spend %3.20f second',t1-t0);
    
    c = lagrange22(A,b,[B;H],[bf;zeros(size(H,1),1)]);
    
    fprintf('\nComputing the L_inf error...'); % adaptive procedure
    err_L2 = posteriori_elliptic(dofs,V,T,d,c,Int_UU2_cell,Int_UV2_cell,Int_UW2_cell,Int_VV2_cell,Int_VW2_cell,Int_WW2_cell);
%     err_L2 = smooth_error(V,posV,T,posT,err_L2);
    error = sum(err_L2);
    fprintf('\nthe total error is %e. matrix size is: %d\n',error,size(A,1));
    if error >= max_err
        flag = 1; % need to comput once more
        err_min = min(err_L2);
        err_max = max(err_L2);
        refine_flag_T = (err_L2 > err_min + 0.6*(err_max - err_min));
        d = d + refine_flag_T;
        count = count +1;
    end
end
toc;  %  show_tri;

%plot the mesh and magnify it.
% figure('Position',[1 250 800 300]);
ZERO=zeros(size(V,1),1); trisurf(T,V(:,1),V(:,2),ZERO);
axis([min(V(:,1)) max(V(:,1)) min(V(:,2)) max(V(:,2))]);
view(0,90);

toc;
%plot error
[Z,Tris,Points] = heval(dofs,V,T,c,d,tri_temp,disp_d);
exact = feval('elliptic_u',Points(:,1),Points(:,2),caseNum);
figure; trisurf(Tris,Points(:,1),Points(:,2),abs(Z-exact));
axis([-1 1 -1 1 0 1e-5]); view(50,30)